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Abstract 

The wave equation for vectors and symmetric tensors in spherical coordinates is studied under 
the divergence-free constraint. We describe a numerical method, based on the spectral decom- 
position of vector/tensor components onto spherical harmonics, that allows for the evolution of 
only those scalar fields which correspond to the divergence-free degrees of freedom of the vec- 
tor/tensor. The full vector/tensor field is recovered at each time- step from these two (in the vector 
case), or three (symmetric tensor case) scalar fields, through the solution of a first-order system of 
ordinary differential equations (ODE) for each spherical harmonic. The correspondence with the 
poloidal-toroidal decomposition is shown for the vector case. Numerical tests are presented using 
an explicit Chebyshev-tau method for the radial coordinate. 
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1. Introduction 

Evolution partial differential equations (PDE) for vector fields under the divergence-free con- 
straint appear in many physical models. Similar problems are to be solved with second-rank 
tensor fields. In most of these equations, if the initial data and boundary conditions satisfy the 
divergence-free condition, then the solution on a given time interval is divergence-free too. But 
from the numerical point of view, things can be more complicated and round-off errors can create 
undesired solutions, which may then trigger growing unphysical modes. Therefore, in the case of 
vector fields, several methods for the numerical solution of such PDEs have been devised, such 
as the constraint transport method [12] or the toroidal-poloidal decomposition [11, 20]. The aim 
of this paper is to present a new method for the case of symmetric tensor fields, which appear in 
general relativity within the so-called 3+1 approach [1], keeping in mind the vector case for which 
the method can be closely related to the toroidal-poloidal approach. We first give motivations for 
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the numerical study of divergence-free vectors and tensors in Sees. 1.1 and 1.2; we briefly intro- 
duce our notations and conventions for spherical coordinates and grid in Sec. 1.3. The case of 
the vector divergence-free evolution is studied in Sec. 2, and the link with the poloidal-toroidal 
decomposition is detailed in Sec. 2.3. We then turn to the symmetric tensor case in Sec. 3 with the 
particular traceless condition in Sec. 3.3. A discussion of the treatment of boundary conditions is 
given in Sec. 4, with the particular point of inner boundary conditions (Sec. 4.3). Finally, some 
numerical experiments are reported in Sec. 5 to support our algorithms and concluding remarks 
are given in Sec. 6. 

1.1. Divergence-free vector fields in relativistic magneto-hydrodynamics 

In classical electrodynamics, the magnetic field is known to be divergence-free since Max- 
well's equations. This result can be extended to general relativistic electrodynamics as well. In 
classical hydrodynamics, the continuity equation can be expressed as d t p + V • (pu) = 0, where p 
is the mass density of the fluid, and u its velocity. Various approximations give rise to divergence- 
free vectors. Incompressible fluids have constant density along flow lines and therefore verify 
that their velocity field u is divergence-free. Water is probably the most common example of an 
incompressible fluid. In an astrophysical context, the incompressible approximation can lead to 
a pretty good approximation of the behavior of compressible fluid provided that the flow's Mach 
number is much smaller than unity. Another useful hydrodynamic approximation is the anelastic 
approximation, which essentially consists in filtering out the sound waves, whose extremely short 
time scale would otherwise force the use of an impractically small time step for numerical pur- 
poses. In general-relativistic magneto-hydrodynamics, the anelastic approximation takes the form 
V • (pTu) = 0, where u is the coordinate fluid velocity, T the Lorentz factor of the fluid, and p its 
rest-mass density. 

Divergence-free vectors have given rise to a large literature in numerical simulations. For 
example, while using an induction equation to numerically evolve a magnetic field, there is no 
guarantee that the divergence of the updated magnetic field is numerically conserved. The most 
common methods to conserve divergences in hyperbolic systems are constrained transport meth- 
ods, projection methods or hyperbolic divergence cleaning methods (see [25] for a review). 

1.2. Divergence-free symmetric tensors in general relativity 

The basic formalism of general relativity uses four-dimensional objects and, in particular, sym- 
metric four- tensors as the metric or the stress-energy tensor. A choice of the gauge, which comes 
naturally to describe the propagation of gravitational waves is the harmonic gauge (e.g. [8]), for 
which the divergence of the four- metric is zero. The 3 + 1 formalism (see [1] for a review) is an 
approach to general relativity introducing a slicing of the four-dimensional spacetime by three- 
dimensional spacelike surfaces, which have a Riemannian induced three-metric. With this formal- 
ism, the four-dimensional tensors of general relativity are projected onto these three- surfaces as 
three-dimensional tensors. Consequently, the choice of the gauge on the three-surface is a major 
issue for the computation of the solutions of Einstein's equations. 

The divergence-free condition on the conformal three-metric has already been put forward 

by Dirac [9] in Cartesian coordinates, and generalized to any type of coordinates in [4]. This 

conformal three-metric obeys an evolution equation which can be cast into a wave-like propagation 
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equation. Far from any strong source of gravitational field, this evolution equation tends to a tensor 
wave equation, under the gauge constraint. With the choice of the generalized Dirac gauge this 
translates into the system we study in Sec. 3, with the addition of one extra constraint: the fact that 
the determinant of the conformal metric must be one (Eq. (167) of [4]). 

The choice of spherical coordinates and components comes naturally with the study of isolated 
spheroidal objects as relativistic stars or black holes. Moreover, boundary conditions for the metric 
or for the hydrodynamics equations can be better expressed and implemented using tensor or 
vector components in the spherical basis. The numerical simulations of astrophysically relevant 
objects in general relativity must therefore be able to deal with the evolution of divergence-free 
symmetric tensors, in spherical coordinates and components. A particular care must be given to 
the fulfillment of the divergence-free condition, since this additional constraint sets the spatial 
gauge on the spacetime. 

1.3. Spherical components and coordinates 

In the following, unless specified, all the vector and tensor fields shall be functions of the four 
spacetime coordinates \(t, r, 9, <p) and h(t, r, 9, (p), where (r, 9, <p) are the polar spherical coordi- 
nates. The associated spherical orthonormal basis is defined as: 

e = 1 ee = -— e = — !— — (1) 

r dr' r 88' v rsin9d(p 

The vector and symmetric tensor fields shall be described by their contravariant components 
[V, V 6 , V*} and [h rr , h rd , h r(p , h 60 , h Gtfi , h^}, using this spherical basis: 

V= J] h = Z Z ^' e i® e i- ( 2 ) 

i=r,9,<p i=r,9,<p j=r,8,tp 

The scalar Laplace operator acting on a field cf>(r, 9, cp) is written: 

d 2 <p 2d(j) 1 
= -t-j + --j- + -A 9ip( f>, (3) 
or 1 r or r l 

where A 6(p is the angular part of the Laplace operator, containing only derivatives with respect to 9 
or cp: 

d 2 (p cos9d<p 1 d 2 <p 
88 1 sin 8 89 sin 9 dip 1 
We now introduce scalar spherical harmonics, defined on the sphere as (see Sec. 18.11 of [2] 
for more details) 

W > 0, Vm, < m < £, Y™(9, <p) = e im(p P'}\cos 9), (5) 
where PJT is the associated Legendre function. For negative m, spherical harmonics are defined 

Vm, -i < m < 0, Y?(0, <p) = (-l) m e im(p P^(cos 9). (6) 

Their two main properties used in this study are that they form a complete basis for the devel- 
opment of regular scalar functions on the sphere, and that they are eigenfunctions of the angular 
Laplace operator: 

W, m), A e(fi Y? = -t{t + 1)17. (7) 
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2. Vector case 



We look for the solution of the following initial-boundary value problem of unknown vector 
V, inside a sphere of (constant) radius R, thus V(0, (p): 

W > 0, Vr < fl, — = AV, (8) 

W > 0, Vr < R, V • V = 0, (9) 
Vr<R, \(0,r,9,<p) = v (r,9,<p), 
d\ 
dt 

W>0, \(t,R,9,<p) = Mt,9,<p). (10) 



Vr<R, 



= w (r, 9, (p), 

f=0 



Vo, Wo and bo are given regular functions for initial data and boundary conditions, respectively. A 
is the vector Laplace operator, which in spherical coordinates and in the contravariant representa- 
tion (2) using the orthonormal basis (1) reads: 

d 2 V 4dV 2V 1 2 
or z r or r L r L r 

(AV) 6 = — + -— + - \Ag v V e + 2— - —j- - 2-^ — 

dr 2 r dr r 2 \ 86 S in 2 9 sin 2 9 dcp 

,*™w d2yv 2dyv 1 (l 2 dyr „cos6dV e V v 

(AYf = —— + -— + — [Ag^ + — - — + 2- 



dr 2 r * r 2 \ sin 9 d<p sin 2 d<£ sin 2 9 

with the divergence 

@ = V-V = + + -I + + . (12) 

or r r \ 89 tan 9 sin 9 dip J 

One can remark that a necessary condition for this system to be well-posed is that 

V • v = V • wo = 0. (13) 

In addition, the boundary setting at r = R is actually overdetermined: the three conditions are not 
independent because of the divergence constraint. This aspect of the problem will be developed in 
more details in Sec. 4.1. 

In the rest of this Section, we devise a method to verify both equations (8) and (9). This 
technique is similar to that presented in [3] with the difference that we motivate it by the use 
of vector spherical harmonics, and can easily be related to the poloidal-toroidal decomposition 
method, as discussed in Sec. 2.3. 

2.1. Decomposition on vector spherical harmonics 

The first step is to decompose the angular dependence of the vector field V onto a basis of pure 
spin vector harmonics (see [24] for a review): 

V(f, r, 9,<p) = J] (E £m (t, r)Yl + B em (t, r)Y B tm + R e "\t, r)Y* m ) , (14) 

e,m 
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defined from the scalar spherical harmonics as 



W>0, V-€<m<e, Yf m = rVyf, (15) 
W>0, M-£<m<£, Y* m = e r xYf m , (16) 
W > 0, V - t < m < £, Y R Cm = Y™ e r ; (17) 

where V is the gradient in the orthonormal basis (1). Note that both Y E tm and Y B tm are purely trans- 
verse, whereas Y R [m is purely radial. From this decomposition, we define the pure spin components 
of V by summing all the multipoles with scalar spherical harmonics (5): 



V(t,r,9,ip) = J^E^Y?, (18) 
W(t,r,9,<p)= Yb^Y?, (19) 



e,m 

the last one being the usual r-component 

Y J R tm Y? = V r . (20) 

l,m 

The advantages of these pure spin components are first, that by construction they can be expanded 
onto the scalar spherical harmonic basis, and second, that angular derivatives appearing in all 
equations considered transform into the angular Laplace operator (7). 

To be more explicit, (V, V M ) can be related to the vector spherical components by (see also [4]): 

V = ^ - -i-f , pi) 
oO smO dip 

yw _ 1 dV | dV» 



sin 6 dip 86 



and inversely 



dV 6 V 9 1 dV v 
d6 tan 6 sin 6 dip ' 
dV v 1 dV e 



\V = — + — + — — (22) 



A ^ = —- + — (23) 
d6 tan 8 sin 9 dip 

Let us here point out that the angular Laplace operator A 0(p is diagonal with respect to the functional 
basis of spherical harmonics and, therefore, the above relations can directly be used to obtain V v 
and V. 

Thus, if the fields are defined on the whole sphere 6 e [0,n], tp £ [0,2^), it is possible to 
transform the usual components \ V 6 , V v ) to the pure spin ones (V, V^) by this one-to-one trans- 
formation, up to a constant {€ = part) for V and V M . Since this constant is not relevant, it 
shall be set to zero and disregarded in the following. Therefore, a vector field shall be represented 

equivalently by its usual spherical components or by (V, V, V 11 ). 
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2.2. Divergence-free degrees of freedom 

From the vector spherical harmonic decomposition, we now compute two scalar fields that 
represent the divergence-free degrees of freedom of a vector. We start from the divergence of a 
general vector W, expressed in terms of pure spin components: 

dW W 1 
= — + 2— + -A e<fi W; (24) 
or r r 

where W has been computed for the vector W from Eq. (22). This shows that the divergence of 
W does not depend on the pure spin component W M . On the other hand, it is well-known that any 
sufficiently smooth and rapidly decaying vector field W can be (uniquely on IR 3 ) decomposed as a 
sum of a gradient and a divergence-free part (Helmholtz's theorem) 

W = V0 + Do, (25) 

with V • D = 0. From the formula (23), one can check that the component W M only depends on 
D . Next, taking the curl of W and, in particular, combining the 9- and <p- components of this curl, 
one has that d r W + ^ - ^ has the same property of being invariant under the addition of any 
gradient field to W, thus depends only on D . Therefore, we define the potential 

dW W W 

A = — — + • (26) 

or r r 

As a consequence, we have that 

D = <=> W M = and A = 0. (27) 

We have thus identified two scalar degrees of freedom for a divergence-free vector field, which 
can be easily related to the well-known poloidal-toroidal decomposition (Sec. 2.3), but have the 
advantage of being generalizable to the symmetric tensor case. 

We now write the wave equation (8) in terms of V and A (computed from V and V). It is 
first interesting to examine the pure spin components of the vector Laplace operator (11): 

(AV)" = AV + 2^-, (28) 
r L 

(AV)" = AV M ; (29) 
one sees that the equation for V decouples from the system, therefore Eq. (8) implies that 

d 2 Vf 



= AV 11 . (30) 

or 
■quation f 

original wave equation (8), we obtain 



8t 2 

Forming then from (11) and (28) an equation for the potential A, which is a consequence of the 



d 2 A 

W = AA. (31) 

We are left with two scalar wave equations, (30) and (31), for the divergence-free part of the 

vector field V. The recovery of the full vector field shall be discussed in Sec. 2.4; the treatment of 

boundary conditions shall be presented in Sec. 4.1. 
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2.3. Link with poloidal-toroidal decomposition 

According to the classical poloidal-toroidal decomposition, a divergence-free vector field F 
can be considered to be generated by two scalar potentials <D and W, via 

F = V x OFk) + V x V x (Ok) (32) 

Here, k is a unit vector, called the pilot vector, which is chosen according to the geometry of the 
problem considered. In [6, 7], k is chosen to be e z in cylindrical coordinates. One can also find the 
decomposition F = V x (A(r, 9)e v ) + B(r, 9)e v when considering axisymmetric solenoidal fields (see 
for example [17]). The latter representation makes A appear clearly as a poloidal component, and 
B as a toroidal component. In order to link the general poloidal-toroidal formalism to our previous 
potentials, we chose k = e r in spherical coordinates (sometimes called the Mie decomposition, see 
[10] ). Then, one can show that 

F = ~A eifi Q> e r + - (—d^ + d e dM e e + - l-d B V + —d v d r Q>\ e v (33) 
r 2 r \ sin 9 J r\ sin 9 J 

Hence, we can identify the former pure spin components F 11 and F 11 through 

F n = -5 r O 
r 

F^ = 

r 

Therefore, the potential A is linked to the potential O via 




which gives us a compatibility condition 

A 0V A = -A(rF r ) (35) 

The latter equation expresses that d r {r 2 &) = for the original vector. Since our vector is a regular 
function of coordinates, it expresses that = 0. 
One can also show the following relations 

e r • V x F = -A 6ip F M 
r 

e r • V x V x F = -AftpA 
r 

2.4. Integration scheme 

We defer to Sec. 5.1 the numerical details about the integration procedure, and we sketch here 
the various steps. From the result of Sec. 2.2, the problem (8)-(9) can be transformed into two 
initial-value boundary problems, for the component V (30) and the potential A (31) respectively. 



Initial data can be deduced from v and w , so that V(/ = 0) and dV 1 /dt(t = 0) are the ji- 
components of, respectively, v and w . The same is true for the A potential. The determination of 
boundary conditions from the knowledge of bo shall be discussed in Sec. 4. We therefore assume 
here that we have computed the component V 1 and the potential A, inside the sphere of radius R, 
for a given interval [0, T], and we show how to recover the whole vector V. 

The pure spin components (V, V) of the vector V are obtained by solving the system of 
PDEs composed by the definition of the potential A (26), together with the divergence-free condi- 
tion (24). From their definitions (18)-(20), it is clear that the angular parts of both V and V can 
be decomposed onto the basis of scalar spherical harmonics, and therefore A as well: 



A(t,r,9,<p)=YA tm (t,r)Y?(e,<p). 



(36) 



We are left with the following set of systems of ordinary differential equations in the r-coordinate: 



( dR 



W > 0, Vm - i < m < I, - 



Cm 



+ 2 ^_f(£+l> £ ,„ = 

dr r r 

— — + = A Cm 

or r r 



(37) 



The potential A being given, the pure spin components V and V are obtained from this system, 
with the boundary conditions discussed in Sec. 4.1. The //-component is already known too, so 
it is possible to compute the spherical components of V V? e [0, T], from Eqs. (21). Note that 
all angular derivatives present in this system (37) are only in the form of the angular Laplace 
operator A 6(p (4). It must also be emphasized that the divergence-free condition is not enforced 
in terms of spherical components (Eq. (12)), but in terms of pure spin components. Thus, if the 
value of the divergence is numerically checked, it shall be higher than machine precision, because 
of the numerical derivatives one must compute to pass from pure spin to spherical components 
(Eqs. (21)). 

The properties of the system (37) are easy to study. Substituting R tm in the first line by its 
expression as a function of E tm and A tm (obtained from the second line), one gets a simple Poisson 
equation: 

A(rE em ) = r + 2A tm . (38) 

\ I dr 

The discussion about boundary conditions, homogeneous solutions and regularity for r = and 

r —* co are immediately deduced from those of the Poisson equation (see e.g. [13]). 

In the case where a source S is present on the right-hand side of the problem (8), the method of 

imposing V-V = can be generalized by adding sources to Eqs. (30)-(31), which are deduced from 

S. Indeed, it is easy to show that the source for the equation for V 1 is the pure spin /i-component 

of S and the source for the equation for A is the equivalent potential computed from S pure spin 

components, using formula (26). Note that an integrability condition for this problem is that the 

source be divergence-free too. Therefore, for a well-posed problem, any gradient term present in 

S can be considered as spurious and is naturally removed by this method, since the /i-component 

and the A potential are both insensitive to the gradient parts. 
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3. Symmetric tensor case 



Similarly to the vector case studied in Sec. 2, we look here for the solution of an initial- 
boundary value problem of unknown symmetric tensor h, inside a sphere of radius R. As explained 
in Sec. 1 .3, the symmetric tensor h shall be represented by its contravariant components h lj (= h jl ), 
where the indices run from l(r) to 3((f); moreover, we suppose that all components of h decay to 
zero at least as fast as 1 /r, as r — > oo. We shall also use the Einstein summation convention over 
repeated indices. 



Thus the problem is written, V(#, ip): 

W > 0, Vr < R, 

W > 0, Vr < R, 
Vr <R, 

Vr <R, 

V? > 0, 



d 2 tij 



= Ati j , 
= 0, 



(39) 
(40) 



h ij (0,r,6,<p) = a i j(r,6, i p), 
dh ij 

t=0 



dt 



h ii (t,R,6,<p)=$(t,6,<p). 



(41) 



The tensors a^y^ and 0^ are given regular functions for initial data and boundary conditions, 
respectively. The full expression of the tensor Laplace operator in spherical coordinates and in the 
orthonormal spherical basis (1) is given by Eqs. (123)-(128) of [4] and shall not be recalled here. 
We point out again that the boundary setting at r = R is overdetermined: this is discussed in more 
detail in Sec. 4.2. 

We introduce the vector H, defined as the divergence of h lj and given in the spherical con- 
travariant components (2) by: 



dh rr 


2h rr 




+ 


dr 


r 


dh r0 






+ 


dr 


r 


dh r * 






+ 


dr 


r 



1 idh 



+ - 



r9 



+ 



1 dh r >? 



r\ 89 sin# dip 
1 (dh 00 1 3h e * 



- h 00 - m + 



h 



r9 



tan6> 



+ - 

r 



+ - 

r 



36 



+ 



+ 



+ 



sin 6 dip tan 
1 dh^ Ih * 



— (h 00 - h™)) , (42) 



+ 



86 smd dip tan# 



= 0. 



We now detail, in the rest of this Section, a method to verify both evolution equation (39) and 
the divergence-free constraint (40). 

3.1. Decomposition on tensor spherical harmonics 

As in the vector case (Sec. 2.1), we start by decomposing the angular dependence of the tensor 
field h lj onto pure spin tensor harmonics, introduced by [21] and [27] (we again use the notations 
of [24]): 



h(t, r, 6, <p) 
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(43) 



h rr (t,r,e,<p) 


\ 1 j tm ym 

l,m 


h T (t,r,9,<p) 


_ Zj 1 i ' 

f,m 

\ 1 T7fm v'" 

= Zj 1 

On 

\ 1 pfm v»i 

~ Zj 1 *' 

_ pta ym 

~ 2 ^ > 


h*(t,r,9,<p) 


W(t,r,e,(p) 


h™{t,r,e,<p) 


h X (t,r,6,p) 





where (z^" 1 , T^"", Ef", B l ™, E e 2 m , B { 2 "*} are all functions of only (t, r). Complete definitions and prop- 
erties of this set of tensor harmonics can be found in [24]. Note that these harmonics have been 
devised in order to describe gravitational radiation, far from any source. In that respect, the most 
relevant harmonics are T El and T Bl , since they are transverse and traceless. The pure spin compo- 
nents of the tensor h are defined as: 

" (44) 

(45) 
(46) 
(47) 
(48) 
(49) 

Explicit relations between the last five components and the usual spherical components (2) are 
now given. 

h T = h 66 + h w (50) 
is transverse; and the total trace is simply given by 

h = h rr + h T . (51) 

In the following we shall use either the component h T or the trace. The components h n and h M have 
similar formulas to those of the vector pure spin components, as {^ r '} ._ 1 2 3 can be seen as a vector: 

(52) 

oO sin0 o(p 

h r<P _ _L^! + ®L. 
sin 9 dp <90 ' 

the reverse formula being similar to Eqs. (22) and (23), they are not recalled here. Finally, the last 
two components are obtained by: 

= (hOO-flW) _ £2^V ^ J_d^ _ J_CMW _ d I 1 dh X \ 

2 36 2 tan 6 88 sin 2 (9 <V 50 \ sin dp]' 

ew _ d 2 h* 1 dh x 1 d 2 h x 2 d ( 1 dh^ \ 
d9 2 tan0 89 s in 2 dtp 2 + 50\sin0 dp ]' 
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and the inverse relations are given by: 



. / . x lW d 2 P 3 dP 

, / , lX d 2 h e * 3 dh e * 
A„,(A„, + 2)^ = _ + __ 



1 d 2 P 



sin 



2 edip 2 



2 d dh e * h 9 * 

IP + + I , 

sin#d<£ \ 86 tan#/ 



sin 



<V 



+ 



sin 9d(f\d9 tan 



(54) 
(55) 



Here as for the vector case, the h n and h M components do not contain any relevant I = term, 
whereas and h x contain neither £ = 0, nor £ = 1 terms, as expected for transverse traceless 
parts of the tensor h. We shall use any set of components of the tensor h: either the usual ones 
\h lj \, using the spherical basis, or the pure spin ones \h rr , h T (or h), hP, \t, ti™, h x }. 

3.2. Divergence-free degrees of freedom 

The vector H defined as the divergence of h in Eq. (42) can be expanded in terms of vector 
pure spin components, which are then written as functions of the tensor pure spin components of 
h (we use the trace h instead of h T ): 



H r = 



dh n 



+ + - i^apff - h) , 

r r v ' 

dh" 3h" \1, v,. 



+ 



h-h r 



W = A H 



dW 3hf 1/ A n s, x 

"7T + + -[Ag v + 2)h X 

or r r v ' 



(56) 
(57) 

(58) 



A possible generalization of the Helmholtz theorem to the symmetric tensor case is that, for any 
sufficiently smooth and rapidly decaying symmetric tensor field T, one can find a unique (on M 3 ) 
decomposition of the form 

T ij = y L J + yj L i + fc y (59) 



with Vjh'J = 0. With these definitions, VjT' J = 



L l = which means that, from the 



six scalar degrees of freedom of the symmetric tensor T' J , the three longitudinal ones can be 
represented by the three components of the vector L. Therefore, the divergence-free symmetric 
tensor h'^ has only three scalar degrees of freedom that we exhibit hereafter. 
One can check that the three scalar potentials defined by 



& = 
3 



dr 

~dr~ 



1 



2r 



TV 



-—l^T w - — + 



J1 J _ jr 



c _dT dT rr + T 
dr dr r 



r 

3T rr 



Ar 



-2A fl 



IfffW T 

+ 



\ dr 



satisfy the property 



JI = S = C = 
11 



h = 0, 



(60) 
(61) 
(62) 

(63) 



and represent the three divergence-free scalar degrees of freedom of a symmetric tensor. 

In order to write the wave equation (39) in terms of these potentials, we first express the pure 
spin components of the tensor Laplace operator acting on a general symmetric tensor h: 

6h rr 4 2h 

(Ah)"" = Ah rr - — - -A g ,h" + — (64) 

AI „ 2dW 2W 2(dh" 3H> , A ^h™ 1, 3 , \ 
(Ah)" = Ah" + -— + — -- — + — + (A ev +2)— + -h- -h rr \, (65) 
r or r L r \ or r v ' r 2r 2r j 

(Ahf = Aff + -— + — — + + (A dlf + 2) — \, (66) 

r or r l r \ or r v ' r ) 

(Ah) w = A^ + ^ + ^, (67) 



y y 2h x 2W 
(Ah)* = Ah x + — + —, (68) 

trace of Ah = Ah. (69) 

The term between parentheses in Eq. (66) is exactly zero in the case of a divergence-free tensor, as 
it represents the //-component of the vector H (58). The similar term in Eq. (65) reduces to -h rr /r, 
when using W = with Eq. (57). We can now write evolution equations, implied by the original 
tensor wave equation (39): 

w = AM, (70) 

The situation is therefore slightly more complicated than in the vector case with Eqs. (30)- 
(31). Indeed, the two potentials !B and C are coupled, but it is possible to define new potentials 
satisfying decoupled wave-like evolution equations. We first write the scalar spherical harmonic 
decomposition of Jl, S and C: 



m, r, e,<p)=J] M {m (t, r)Y?{6, <p), 

e,m 

8(t,r,0,<p) = Yj&^rWfid,?), 

e,m 

C(t,r,0,<p)= YjC^it^ieM- 
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Then, we define new potentials S and C as: 



S(t, r, 0,tp)=Yi ( 2S >> r) + Y?(9, <p), 

C(t, r, 9,<p)=J] r ) ~ ^ m (t, rj) Y?(9, <p). 

The Eqs. (71)-(72) are transformed into: 



(73) 
(74) 





d 2 S 




dt 2 




d 2 C 




dt 2 


with, for any scalar field f(r, 6, <p) = 


Z(o«) f tm (r 


d 2 f 


101 + L 

r dr r 2 ' 


d 2 f 


28f 1 
r dr r 2 ' 



= AS, 
= AC; 



(75) 
(76) 



tm 



Im 



(77) 
(78) 



These two operators are very similar to the usual Laplace operator, but in the angular part Ag v , 
they contain a shift of, respectively -1 and +1 in the multipolar number I, for A and A. We thus 
have obtained three evolution wave-like equations (70), (75) and (76) for the three scalar degrees 
of freedom of a divergence-free symmetric tensor. 

3.3. Traceless case 

As presented in Sec. 1.2, some evolution problems of symmetric tensors in general relativity 
can have another constraint, in addition to the divergence-free condition already studied (40). 
This is the condition of determinant one for the conformal metric which turns into an algebraic 
condition (Eq. (169) of [4]), and is enforced by iteratively solving a Poisson equation with the 
trace of the tensor as a source, as described in Sec. V.D of [4]. Therefore, in the following the 
trace of the unknown tensor h is assumed to be known. 

The fact that the trace h (51) of a divergence-free symmetric tensor is fixed reduces a priori 
the number of scalar degrees of freedom to two. For instance, we here show that if the trace is 
given, the scalar potentials S and C are linked. We take the partial derivative with respect to r of 
the definition of C (62) and B (61) to obtain: 



dC 2C 2A IdS 3S _ C , , 
dr r e<p \ dr r Arj 



(79) 



Therefore, if h and C are given, it is possible to integrate this relation with respect to the r- 

coordinate to obtain S (which we have assumed to converge to as r —> oo). Because of the 

definitions (73)-(74), S and C are also linked together if the trace is given. 
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We shall assume in the following that this trace is zero. All the equations presented hereafter 
can easily be generalized to the non-zero (given) trace case, taking the general form of the equa- 
tions of Sec. 3.2. We shall therefore use only two scalar potentials, namely J?l and S to describe a 
general traceless divergence-free symmetric tensor. 

3.4. Integration scheme 

Similarly to what has been done in the beginning of this section, we consider the homogeneous 
wave equation for a symmetric tensor (39), under the constraints that the tensor be divergence- 
free (40) and traceless (h = 0). We have seen in Sec. 3.3 that it was necessary to solve for at least 
the two wave-like evolution equations (70) and (75). We describe now how to obtain the whole 
tensor, once !Mt, r, 6, ip) and S(t, r, 9, ip) are known. 

In order to obtain first the six pure spin components (actually, their spherical harmonic decom- 
positions (44)-(49)) of h at any time t, we use the following six equations: the traceless condition, 
the three divergence-free conditions and the definitions of J{ and B. They represent two systems of 
coupled differential equations in the r-coordinate, that we express in terms of the tensor spherical 
harmonic components (43). The first one comes from the definition of Jl (60) and the W = 
condition (58); it couples the [i- and the /V-components of h: 

dB* m B\ m 

— = rf m , (80) 

or r 

dB\ m 3B e ' n 2-t(t+X)Bi m 

^ + — + ^-=0. (81) 

or r r 

This system has two unknown functions B l ™ and B%", whereas Jl em is obtained from the time 
evolution of 3\(t, r, 6, tp). 

The second one comes from the definition of & (73) and the two W = W = conditions (56)- 
(57); it couples the rr-, rj- and "W-components: 

dE tm E tm 2E\ m 1 dL e n m £ + 4Li m 

(I + 2)— 2- + 1(1 + 2)-2 !- - — — - - — = S {m , (82) 

or r r 2(1 + 1) or i + 1 2r 

dLi m 3Li m 1(1 + \)E\ m 

-^ + — — = 0, (83) 

or r r 

dE lm 3E em L tm 2 - £(£ + l)E £ ™ 

-T- + — ~ -T- + — = 0. (84) 

or r 2r r 

Here, the unknowns are L e m , Ef" and E e 2 m and & m is known from the evolution of S(t, r, 6, <p). 

When looking at a more general setting, the trace h appears only in the second system. If we 
combine Eq. (80) with Eq. (81), we obtain a Poisson equation for the unknown rh x , with J\ and its 
radial derivative as a source. As for the vector case, this system can be solved using, for example, 
the spectral scalar Poisson solver described in [13], and one obtains the pure spin components h? 
and h x . 

Such an argument cannot be used for the second system, but a search for homogeneous solu- 
tions gives that, for a given I, the simple powers of r. 

1 1 and -L (85) 



r {+3 r e+\ 
r 14 ' 



represent a basis of the kernel of the system (82)-(84). With this information, one can devise a 
simple spectral method to solve this system (see Sec. 5.1) and obtain the pure spin components 
h rr ,h! 1 and h^. With the traceless condition, one can also recover h T from h rr , and finally use 
Eqs. (52)-(53) to get the spherical components of h. 



4. Boundary conditions 

4. 1 . Vector system 

We discuss here the spatial boundary conditions to be used during our procedure, so that we 
recover the unknown vector field at any time-step. The source of the vector wave equation is put 
to zero for the sake of clarity; but the reasoning would be exactly the same in the general case. 

As pointed out in Sec. 2.4, the recovery of the vector field at each time-step will require two 
different operations: first, we use the two scalar wave equations (31) and (30) to recover A and 
V. Two boundary conditions, set at the outer sphere (the boundary of our computation domain), 
will then be needed for these quantities. The second step will consist of the inversion of the 
differential system (37), to obtain the pure spin components V and V. This system is, in terms 
of the structure of the space of homogeneous solutions, mathematically equivalent to a Poisson 
problem (see Eq. (38)); its inversion will then also require an additional boundary condition. 

From the setting of our problem presented at the beginning of Sec. 2, we can impose Dirichlet 
boundary conditions for the 3 pure spin components on the outer sphere. The condition on V 
enables us to recover the value of the entire field on our computational domain, through the direct 
resolution of (30). Once we obtain the value of the field A on our domain, we can use a condition 
on either V or V to invert the system (37), and retrieve the additional spin components. 

There remains the necessity of imposing a boundary condition on A to solve Eq. (31). This 
cannot be done using condition at r = R in (10) and the definition (26), because ^ must be 
specified. To overcome this difficulty, we exhibit here algebraic relations that link the value of A 
at the boundary and time derivatives of the pure spin components. These will be compatibility 
conditions, derived only from the structure of our problem. We express radial derivatives of equa- 
tions (24) and (26), respectively, to obtain, using relations (11) and (28), the following identities 
(see also Eq. (35)): 

1 d 2 V r 

A A = "Mr- (86) 
M + A _ 

dr r dt 2 

Those equations are derived using only the fact that our vector field satisfies the wave equation 
and is divergence-free. From the knowledge of the vector field at the boundary, we can impose 
either of these two relations as boundary conditions for A; the first being of Dirichlet type for each 
spherical harmonic of A, the second of Robin type. This way we are able to solve equation (31), 
and complete our resolution scheme. 

Let us finally note that our boundary problem is, as one could guess, actually overdetermined: 

there is no need to know the value of the entire vector field on the outer sphere. It can be easily 

seen that, if one only has access to the boundary values of V M and V, or V M and V, the boundary 
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conditions for all equations can be provided. This also gives us insight about what would happen 
if we set up a numerical problem in which spatial boundary conditions are not consistent with a 
solution of Eqs. (8, 9); this could occur for example because of numerical rounding errors or simply 
a physical boundary prescription which is not compatible with a divergence-free vector field. Our 
method will then still provide a solution that is divergence-free and which satisfies Eqs. (8, 9); 
however only the boundary conditions that are directly enforced will be satisfied. For example, 
if we choose in our scheme to enforce boundary conditions on V M and V n , the outer boundary 
conditions that are satisfied at each time-step are actually of the form (we keep the notation of 
(10)): 

W>0, V»(t,R,0,<p) = lf Q (t,9,<p), 

y(t,R,e,<p) = bp,e,<p), 

gVg^ftjg) + 2 Q (p) = _I A b v M(pl (88) 

or r r 

The last condition is directly derived from the vanishing of the divergence (Eq. (24)) at the 
boundary. Let us note that we do not even impose a Dirichlet condition on V as was originally 
intended. We may then not satisfy all the boundary conditions we wished to prescribe at first. This 
may also depend on the boundary value we choose to use for the inversion of the system (37). 

We do not treat alternative cases for the boundary problem (for which the knowledge of the 
vector field on the outer sphere could be substituted by, for example, the knowledge of its first ra- 
dial derivative); but a similar approach would also provide expressions for the boundary conditions 
of all the equations tackled in our scheme. 

4.2. Tensor system 

The tensor problem presents itself in a similar way to the vector case, only with a few additional 
difficulties. As seen in Sec. 3.4, we can separate the problem into two parts; the first consists in 
retrieving the field tft from Eq (70), and then get the spin components hf and h x . In a similar way, 
we compute the value of !B from Eq. (75), so that we obtain the fields h rr , and W 1 from the 
inversion of the system (82, 83, 84) . The field h T is deduced from the traceless hypothesis. The 
tensor field is then entirely determined. 

As in the vector case, the solution of wave equations for 3\ and S requires one boundary 
condition for each equation. The elliptic system (80, 81) is also quite similar to that for the vector 
case, and its space of homogeneous solutions is also equivalent to that of a single Poisson equation. 
One boundary condition is also required; it will be chosen as a Dirichlet condition on either h M or 
h x , according to the setting of our problem (41). 

For the elliptic system (82, 83, 84), the homogeneous solutions have been characterized in 
Sec. 3.4. The only basis vector of the kernel of solutions that is regular in our computation domain 
is, for any € > 2, the solution r [ ~ 2 . The other two vectors of the kernel basis are not regular at 
the origin of spherical coordinates. This means, from a basic point of view, that one boundary 
condition will be sufficient at the outer sphere. It will be provided, again according to our problem 
setting, as a Dirichlet condition on any of the fields h rr , hP or 
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The last boundary problem concerns the fields J?l and !B. They will be handled the same way as 
in the vector case. We take the radial derivatives of the equations (58) and (60), using the elliptic 
equations (66) and (68), to obtain the following compatibility conditions: 



d 2 ^ 



dt 2 

^ + 2— - 

dr r dt 2 



(89) 
(90) 



These are again derived using only the divergence-free property of the vector field as well as 
the verification of the main wave equation. Using the known value of, respectively, U 1 and h x at 
the outer boundary, we obtain either a Dirichlet boundary condition for each spherical harmonic 
from the first relation, or a Robin condition with the second one. Again those identities have been 
obtained only from the equations of our problem and the definitions of the variables we use. 

Taking the same path for the second part of the problem, we express radial derivatives of 
Eqs. (56), (57), (61) and (62) to obtain respectively, and for each spherical harmonic, the follow- 
ing relations: 



d 2 L e m 
dt 2 

d 2 E\ m 
dt 2 

d 2 E e 2 m 
dt 2 

d\lJ™ + T e Q m ) 



1 



(2t + l)r 
1 



(2t+\)r 



{t + l)(t- l)S em + 



t + 2 



-C 



Cm 



1 



2^+1 
1 



2^+1 



(£+l)dS {m \dC im (£ + + 2) B £m £-3C e " 
2 d7~ ~ A~r 2 r 4 7 

{€+ l)(£ + 2)dC {m 



(91) 
(92) 
(93) 



8r 



- tit + l)(t + 2)—- + tit + \){t - 1) 



dr 



1 Q£m 

+ -it+l)[tit-3) + t + 4] 

2 r 



(94) 



When expressing the vanishing of the trace, the last equation can be transformed into: 

1 



d 2 E l 2 m 
dt 2 



2t{t+l){2t+ 1) 



dr 



dr 



- t{t+l){t-3)- 



(95) 



Although those equations involve both the fields & and C, one can easily see that combining 
them can lead to conditions on the field S only. For example, the combination of (91) and (92) 
provides, for each index t: 



(t+W-l) 



d 2 Ll m d 2 E\ m 
+(7+1) 



dt 2 



dt 2 



(96) 
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which is interpreted as a Dirichlet boundary condition for S. Robin boundary conditions can be 
obtained from the combination of Eqs. (93), (94), and either (91) or (92). The tensor bound- 
ary problem is then entirely solved; tests for some of the boundary conditions derived here are 
presented in Sec. 5. Let us note again that this problem is overdetermined: concerning the first 
system, the knowledge of a Dirichlet condition on either only hf, or only h x suffices to provide 
boundary conditions for and the system (80, 81). For the part of the algorithm related to S, we 
easily see that Dirichlet conditions for any two of the spin components h rr , h 11 and are sufficient 
to solve the boundary problem. 

We finally point out that, in the same fashion as in the vector case, if the value imposed 
as a Dirichlet condition for the tensor at the outer boundary (Eq. (41)) is not consistent with the 
system, the boundary conditions actually imposed on our scheme will be slightly different: only 
the Dirichlet conditions for the pure spin components that are explicitly enforced will be satisfied. 
Other boundary values will only express the coherence with respect to the fact that the solution 
is indeed divergence free. As done in Sec.4.1, it is possible to express other boundary conditions 
enforced in practice by using the expression for the tensor divergence H as a function of the pure 
spin components. 

4.3. Working in a shell: inner boundary conditions 

We say a few words here about the resolution of the tensorial problem when our computation 
domain is no longer an entire sphere, but is instead bounded on the interior at a finite coordinate 
radius r = R in > 0. We add in our setting the condition that, V(0, tp): 

w>o, h i \t,R in ,e,<p) = ^(t,e,<p). 

Physical information is then also provided at the internal boundary (this is, again, an overdeter- 
mined set of boundary conditions). This new geometry will imply the need for two inner boundary 
conditions to be imposed for the wave equations on 3i and in S. These are easily found using the 
results of the last section and the knowledge of Dirichlet boundary conditions on the inner and 
outer sphere for all components. The system (80, 81) also needs an additional (inner) boundary 
condition, imposed on either W or h x . There is, however, a slight subtlety concerning the triple 
system (82, 83, 84). As seen in Sec. 3.4, the kernel of solutions to this system is of dimension 3, 
and since our computational domain no longer includes r = 0, all 3 basis vectors of this kernel 
are regular in our domain. This means that 3 boundary conditions have to be imposed overall 
for inverting this system (in contrast with the sphere case, where we only imposed one). Those 
three conditions are imposed here on either h rr , W or on each limit of the domain. We have a 
priori the freedom to choose which boundary conditions we want to impose, and where to impose 
them; numerical experimentation would be required to indicate whether or not there are preferable 
choices. 

To conclude this section, we mention also the work of [26] where the authors used the for- 
malism presented in this paper to solve a tensor elliptic equation that is part of a formulation of 
the Einstein equations. The resolution was made on a 3-space excised by a sphere of fixed co- 
ordinate radius, where the tensor equation possessed a weak singularity property (see [15]). The 
boundary condition problem was treated a little bit differently, as all boundary conditions imposed 
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were either emanating from the very structure of the problem, or were not needed at all. This is 
a consequence of the particular behavior of that operator at the boundary; on this setting for the 
domain geometry, one boundary condition was imposed to invert the system in h M and h x , and two 
for the system involving h rr , h n and h™. 



5. Numerical tests 

5.1. Spectral methods in a sphere 

The numerical schemes presented in previous sections have been implemented using a multi- 
domain spectral method in spherical coordinates (see e.g. [2, 16], for general presentations and [14] 
for a more detailed description in the case of numerical relativity). We have used the lorene nu- 
merical library [19], with scalar fields decomposed onto a basis of Chebyshev polynomials, in 
several domains, for the r-coordinate, Fourier series for the ^-coordinate and either Fourier or 
associated Legendre functions for the ^-coordinate (Pf(cos6), see Sec. 1.3). This last option is 
obviously needed by our algorithms, which strongly rely on spherical harmonics decompositions 
and on the angular part of the Laplace operator A 6(p . The other basis of decomposition (Fourier) 
is quite useful for computing angular derivatives d/d6 and operators such as 1/ sin#, appearing 
in e.g. (21) or (52). The coordinate singularity on the z-axis (6 = 0, 7r) is naturally handled by 
the spherical harmonic decomposition basis. We cope with the coordinate singularity at the origin 
(r = 0), using an even/odd radial decomposition basis (only even/odd Chebyshev polynomials), 
depending on the parity of the multipole I (see [5] and Sec. 3.2 of [14]). The complete regularity 
requirement would be that, for each multipole I the radial Taylor expansion of a regular function 
should include only r p with p > i. We have found however that the simpler parity prescription 
described above is in practice sufficient for the study of the wave or Poisson equations performed 
here. 

The wave equations (30)-(31) and (70)-(75) are integrated numerically by writing them as 
first-order systems: 

d 2 <p 



% < 97) 

Tt-^ 

After discretization in the angular coordinates using spherical harmonics, we then use a third-order 
Adams-Bashforth (explicit) time-stepping scheme with a fixed time-step dt and a Chebyshev- 
tau technique in the radial coordinate. The differential systems for the computation of pure spin 
components from the divergence-free degrees of freedom, as system (37) in the vector case, or 
systems (80)-(84) in the tensor case, are solved at every time-step in the Chebyshev coefficient 
space. A tau method is used to match together the solutions across the domains, and to impose the 
boundary conditions at r = R. 

5.2. Vector wave equation 

We consider here the numerical solution of the problem (8)-(10), with v l (r, 6, <p) given by its 
Cartesian components by (with z = rcos(#)): 

vl = -v y = cos(z), (98) 
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Figure 1 : Decay of the errors (difference with theoretical solution and divergence of the numerical solution) for the 
vector wave equation, as a function of the number of radial Chebyshev coefficients N r used in each domain. Other 
settings are R = 6, dt = 0, 00032, N =n,N ip = 4. 

the other component is zero. Thus, the vector v' Q is clearly divergence-free. With appropriate 
boundary conditions, the solution of the problem (8)-(10) is (still in Cartesian components) simple 
to express: 

V\t, r, 9, <p) = -V y (t, r, 9, <p) = cos(0 cos(z), (99) 

the other component being zero. The vector wave equation is solved through the two scalar wave 
equations for the potentials A and the component V M as explained in Sec. 2.4. From Eq. (99), 
we know the values of b l Q (t,6,<p) appearing in Eq. (10) as Dirichlet boundary conditions and we 
can deduce its pure spin components (b^b^bfy. These are used to obtain Dirichlet boundary 
conditions for the evolution equations for A and fi, as described in Sec.4.1 using Eq. (86) for A. 
Finally, the elliptic system (37) is solved with the appropriate Dirichlet boundary condition given 
by the spin component b r (see also Sec. 4.1). 

We use the numerical techniques given in Sec. 5.1, with two domains, and numbers of points 
in each direction given by (N r , Ng, N^. We have integrated the vector wave equation over the time 
interval t e [0, 2n] and looked at the maximum in time of two quantities to estimate the accuracy of 
the solution. First, the difference between the numerical solution and the theoretical one (99), ro- 
tated to spherical basis (1), is computed. Then, the divergence of the numerical solution, expressed 
in the spherical basis is also monitored. Note that, even though all the Cartesian components of V 
do not depend on the azimuthal angle ip, the spherical components do depend on <p and we have 
always used four points in the ^-direction. 

In Fig. 1, we observe as expected an exponential convergence of both the discrepancy between 
the theoretical and numerical solutions (maximum over all grid points and all components) as 
functions of the number of spectral coefficients used in the radial direction N r , all other parameters 
being fixed. The same behavior has been observed when keeping N r fixed and varying N e . Besides, 
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Figure 2: Decay of the errors (difference with theoretical solution and divergence of the numerical solution) for the 
vector wave equation, as a function of the time-step dt. Other settings are R = 6,N r = 11, Ng- \1,N V =A. 



we observe an exponential decay of the divergence of the solution in the second (or outer) domain, 
whereas the divergence of the solution in the first (central) domain remains constant to the radial 
precision. This is due to the matching across domains and imposition of boundary conditions, 
which can be seen as a modification of the solution of the system (37) by the addition of a linear 
combination of homogeneous solutions. These homogeneous solutions of the system (37) are, for 
each multipole t, r e ~ x and 1 // +2 . The latter being singular at r = is not relevant in the central 
domain. The r e ~ l function is a polynomial and is well represented in the first domain, whereas in 
the second domain, we also need to resolve 1// +2 , which is poorly approximated for low values 
oiN r . 

On the other hand, when varying the time-step dt, the difference between the numerical and 
exact solutions decreases as 0(dt 3 ) (see Fig. 2), as expected for a third-order scheme. Another fea- 
ture verified in Fig. 2 is the fact that the divergence of the solution is (almost) independent of the 
time-step, being thus only a function of the spatial resolution. The best accuracy observed in Fig. 1 
is limited by angular resolution and the fact that the divergence is computed using spherical com- 
ponents (Eq. 12), whereas the divergence-free constraint is imposed using pure spin components 
(Eq. 24). Therefore, the computation of derivatives in Eqs. (21) to obtain the spherical components 
introduces additional numerical noise, depending on the angular resolution. 

5.3. Divergence-free and traceless tensor wave equation 

Similarly to Sec. 5.2, we consider here the numerical solution of the problem (39)-(41), with 
a l Q (r, 8, (p) given in the Cartesian basis by (with z = rcos(9)): 

a x Q x = -c% = cos(z), (100) 
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Figure 3: Decay of the errors (difference with theoretical solution and divergence of the numerical solution) for the 
tensor wave equation, as a function of the number of radial Chebyshev coefficients N r used in each domain. Other 
settings are R = 6,dt = 0.00032, N g = 17, = 4. 

all the other components are zero. Thus the tensor dl is clearly symmetric, divergence-free and 
trace-free. With yl = and appropriate boundary conditions, the solution of the problem (39)-(41) 
is (still in Cartesian components) simple to express: 

h xx (t, r, 9, (p) = -h yy {t, r, 9, <p) = cos(0 cos(z), (101) 

all the other components being zero. The tensor wave equation is solved through the two scalar 
wave-like equations for the potentials Jl and S as explained in Sec. 3.4. From Eq. (101), we know 
the values of /ftf(f, 9, <p) appearing in Eq. (41) as Dirichlet boundary conditions and we can deduce 
its pure spin components (j3^ These are used to obtain Dirichlet boundary conditions for 

the evolution equations for 3\ and S, as described in Sec. 4.2 using Eqs. (89) and (96), respec- 
tively. Finally, the elliptic systems (80)-(84) are solved with the appropriate Dirichlet boundary 
conditions given by the spin components of 0^ , namely and We have integrated the tensor 
wave equation following the same procedure as in Sec. 5.2. results are displayed in Figs. 3 and 
4, where we observe as expected an exponential convergence of both the discrepancy between 
the theoretical and numerical solutions, and the divergence of the numerical, as functions of N r . 
When varying the time-step dt, the difference between the numerical and exact solutions decreases 
as 0(dt 3 ), as expected. Here again, the divergence of the solution is (almost) independent of the 

time-step, being thus only a function of the spatial resolution, from the same reasons as in the 
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Figure 4: Decay of the errors (difference with theoretical solution and divergence of the numerical solution) for the 
tensor wave equation, as a function of the time-step dt. Other settings are R-6,N r - 17, Ng = 17, N v = 4. 

vector case. 

6. Concluding remarks 

We have described a new numerical method for solving the wave equation of a rank-two sym- 
metric tensor on a spherical grid, ensuring the divergence-free condition on this tensor. In order 
to describe this method, we have first addressed the vector case, for which we have reformulated 
the poloidal-toroidal decomposition in spherical components. This approach, which relies on a 
decomposition onto vector spherical harmonics was then generalized to the case of a symmetric 
tensor. Through numerical tests of the vector and tensor wave evolution in a sphere using spectral 
explicit time schemes, we have observed that this method was convergent and accurate. In partic- 
ular, the level at which the divergence-free condition is violated is determined only by the spatial 
discretization and does not depend on the time-step, as expected. This method strongly relies on 
the decomposition onto spherical harmonic spectral bases, but is not bound to spectral methods 
for the representation of the radial coordinate. 

The discussion in Sec. 4 gave us the compatibility conditions (86), (89) (96), which are nec- 
essary to obtain boundary conditions for the additional scalar field equations, representing the 
evolution of the divergence-free degrees of freedom of our objects (A, Jl, &). The numerical tests 
performed in this study have dealt only with simple Dirichlet boundary conditions. However, it 
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would be rather straightforward to generalize them to more complex boundary conditions, which 
are needed in realistic simulations of gravitational waves [18, 22, 23]. 

In this respect, an interesting issue would probably be the general well-posed nature of these 
boundary conditions with respect to our scheme, and how the modifications for these conditions 
with this method, sketched in Sec.4. 1 and 4.2, would alter the physical behavior of the solution. 
One could for example think of a Robin-like boundary setting linked to an outer wave-absorbing 
condition (as in [22]), instead of the Dirichlet setting studied here; the fact that boundary condi- 
tions may be only partially verified could have an effect on how this required feature at the bound- 
ary would be described eventually in our scheme. The same type of questions arise in a more 
general case, where the source terms of the equations are non- vanishing: these sources would 
also require well-posedness conditions (i.e. a vanishing divergence for the wave equation). If this 
requirement is not satisfied (because of the iteration procedure or numerical errors), although the 
problem is then mathematically ill-posed, our scheme will still converge: it provides us with a so- 
lution of the wave equation with a source that is basically the divergence-free part of the original 
ill-posed source. The influence of this feature on the general stability and physical relevance of 
the procedure is an open issue. 

Future studies include the simulations of perturbed black hole spacetimes, with the extraction 
of gravitational waves, and the solution of general-relativistic magneto-hydrodynamics in the case 
of a rotating neutron star. 
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